#ifndef __FUNCTIONS_H__
#define __FUNCTIONS_H__

#include "FunctionSet.h"

// Boundary Functions
Real bottom(const Point<2> &p)
{
    return p[1];
}

Real up(const Point<2> &p)
{
    return 1 - p[1];
}

template <int N>
Real left(const Point<N> &p)
{
    return p[0];
}

template <int N>
Real right(const Point<N> &p)
{
    return 1 - p[0];
}

template <int N, int C>
Real constant(const Point<N> &p)
{
    return C;
}

//-------------------- 1D FUNTIONS --------------------
// 1. u(x) = exp(x)
Real f11(const Point<1> &p)
{
    return -exp(p[0]);
}

Real u11(const Point<1> &p)
{
    return exp(p[0]);
}

FunctionSet<1> F11 = {
    left<1>, right<1>,
    constant<1, -1>, u11,
    f11, u11};

// 2. u(x) = x^2

Real u12(const Point<1> &p)
{
    return p[0] * p[0];
}

FunctionSet<1> F12 = {
    left<1>, right<1>,
    constant<1, 0>, constant<1, 2>,
    constant<1, -2>, u12};

// 3. u(x) = sin(\pi x)

Real u13(Point<1> p)
{
    return sin(M_PI * p[0]);
}

Real f13(Point<1> p)
{
    return M_PI * M_PI * sin(M_PI * p[0]);
}

Real u_p(Point<1> p)
{
    return -M_PI;
}

FunctionSet<1> F13 = {
    left<1>, right<1>,
    u_p, u_p,
    f13, u13};

FunctionSet<1> *funcs_1d[3] = {&F11, &F12, &F13};

//-------------------- 2D FUNTIONS --------------------
// 1.u(x,y) = exp(y+sin(x))

double f21(const Point<2> &p)
{
    double x = p[0], y = p[1];
    return (sin(x) - cos(x) * cos(x) - 1) * exp(y + sin(x));
}

double u21(const Point<2> &p)
{
    double x = p[0], y = p[1];
    return exp(y + sin(x));
}

double u21_x(const Point<2> &p)
{
    double x = p[0], y = p[1];
    return cos(x) * u21(p);
}

double u21_y(const Point<2> &p)
{
    return u21(p);
}

FunctionSet<2> F21 = {
    bottom, up, left<2>, right<2>,
    [](Point<2> p)
    { return -u21_y(p); },
    [](Point<2> p)
    { return u21_y(p); },
    [](Point<2> p)
    { return -u21_x(p); },
    [](Point<2> p)
    { return u21_x(p); },
    f21, u21};

// 2.u(x,y) = x^2+y^2

double u22(Point<2> p)
{
    return norm(p) * norm(p);
}

FunctionSet<2> F22 = {
    bottom, up, left<2>, right<2>,
    constant<2, 0>, constant<2, 2>, constant<2, 0>, constant<2, 2>,
    constant<2, -4>, u22};

// 3.u(x,y) = (x^2-x^4)(y^4-y^2)

Real f23(const Point<2> &p)
{
    return 2 * ((1 - 6 * p[0] * p[0]) * (std::pow(p[1], 2) - std::pow(p[1], 4)) + (1 - 6 * p[1] * p[1]) * (std::pow(p[0], 2) - std::pow(p[0], 4)));
}

Real u23(const Point<2> &p)
{
    Real x = p[0], y = p[1];
    return (x * x - std::pow(x, 4)) * (std::pow(y, 4) - y * y);
}

FunctionSet<2> F23{
    bottom, up, left<2>, right<2>,
    constant<2, 0>, constant<2, 0>, constant<2, 0>, constant<2, 0>,
    f23, u23};

FunctionSet<2> *funcs_2d[] = {&F21, &F22, &F23};

#endif // __FUNCTIONS_H__